
**************************************************************************************
**************************************************************************************
*********************************** MAIN TABLES **************************************
**************************************************************************************
**************************************************************************************

* NOTE: All file paths and data names that need to be inserted are identified with <>

clear
set more off
capture log close 

global dta = "<insert path to raw data>" 
global output= "<insert path to outputs>" 

log using "${output}\log\create_tables.log", replace

di "Time $S_DATE $S_TIME"

/* Programs that need to be installed:
rdplot 
regsave 
*/

**************************************************************************************

**************************************************************************************
************************************* TABLE 1 ****************************************
**************************************************************************************
*** The Effect of Disability Insurance on Mortality Rates

local tbl table1

local sample main

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'.dta", replace
restore

*** Generating the regression results

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* Fuzzy RK without covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3)  

  *first-stage results
  local f_b_conv = e(tau_T_cl)
  local f_se_conv = e(se_tau_T_rb)
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)

  *main results - converted to p.p. change per $1000/year
  local b_conv = e(tau_cl) * (1/12) *(1000)* 100
  local se_conv = e(se_tau_cl) * (1/12) * 1000  * 100
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Conv coefficient: " "`b_conv'"
  di "Conv std error: "   "`se_conv'"
  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Conventional 
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(f_conv, `f_b_conv', `f_se_conv', m_conv, `b_conv', `se_conv', bw_conv, `bw', miss, obs_conv, `obs', miss) append

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid /// 
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append


* Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *first-stage results
  local f_b_conv = e(tau_T_cl)
  local f_se_conv = e(se_tau_T_rb)
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)

  *main results - converted to p.p. change per $1000/year
  local b_conv = e(tau_cl) * (1/12) *(1000)* 100
  local se_conv = e(se_tau_cl) * (1/12) * 1000  * 100
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Conv coefficient: " "`b_conv'"
  di "Conv std error: "   "`se_conv'"
  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Conventional 
regsave using "${output}/tables/`tbl'.dta", autoid /// 	
	addvar(f_conv_cov, `f_b_conv', `f_se_conv', m_conv_cov, `b_conv', `se_conv', bw_conv_cov, `bw', miss, obs_conv_cov, `obs', miss) append

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append


*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14 

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen panel = .
replace panel = 1 if bendpoint=="lower"
replace panel = 2 if bendpoint=="fm"
replace panel = 3 if bendpoint=="upper"

gen obsnum = _n
sort panel _id obsnum
drop obsnum

cap gen column = _id
replace column = _id if column==.
replace column = 5   if column==. 
drop _id

save "${output}/tables/`tbl'.dta", replace
}

order var
order panel
order bendpoint

save "${output}/tables/`tbl'.dta", replace


**************************************************************************************
************************************* TABLE 2 ****************************************
**************************************************************************************
*** Table 2 Heterogeneity in the mortality effects at the lower and family bend points

local sample main
local bplist lower fm
foreach bp of local bplist  {

**************************************************************************************
*** Full sample

local tbl table2_`bp'_all

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef = .
gen categ = ""
save "${output}/tables/`tbl'.dta", replace
restore

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

*** Variable descriptions - 0/1 identifiers

* hearings: 0 "DDS"	 		1 "hearings"

* male: 	0 "females"		1 "males"

* black: 	0 "nonblack"	1 "black"

* aged45:	0 "aged<45"		1 "aged 45+"
cap drop aged45
gen aged45 = 1
replace aged45 = 0 if age_filing>=1 & age_filing<45
tab aged45

* yrs0609:	0 "1997-2005"	1 "2006-2009"
cap drop 
gen yrs0609=0
replace yrs0609=1 if startpay_year>=2006 & startpay_year<=2009
tab yrs0609

*** Variable descriptions - other identifiers

* primary disability: individual dummy variables identifying them
* "mental"
* "musc"
* "neoplasms"
* "circulatory"
* "other"

*** RKD with covariates for 0/1 identifiers

local vlist "hearings male black aged45 yrs0609 mental musc neoplasms circulatory other"
foreach v of local vlist  {	

forvalues i=0/1	{
	
preserve
keep if `v'==`i'

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"
	
replace categ = "`v'_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace

restore

}

}
	
*** Mortality in first 2 years

foreach i in 12 34	{
 
use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* time on DI - years 1&2, years 3&4
cap drop d_year12
gen d_year12 = (d_year1 + d_year2) / 2
cap drop d_year34
gen d_year34 = (d_year3 + d_year4) / 2

* Fuzzy RK with covariates 
capture noisily rdrobust d_year`i' ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"

regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year`i'

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

replace categ = "mortality_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}

*** Adding column and group identifiers

use "${output}/tables/`tbl'.dta", clear

cap drop N
cap drop if var=="RD_Estimate"
cap drop if coef==.
cap drop _id

cap gen column=1
replace column=2 if var=="mean" 

gen obsnum = _n
gsort column -obsnum

drop obsnum

order categ
order column

save "${output}/tables/`tbl'.dta", replace

**************************************************************************************
*** DDS only

local tbl table2_`bp'_dds

local cov_controls = "age_filing male black musc mental"

*** Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef = .
gen categ = ""
save "${output}/tables/`tbl'.dta", replace
restore

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* DDS restriction
keep if hearings==0

*** Variable descriptions - 0/1 identifiers

* hearings: 0 "DDS"	 		1 "hearings"

* male: 	0 "females"		1 "males"

* black: 	0 "nonblack"	1 "black"

* aged45:	0 "aged<45"		1 "aged 45+"
cap drop aged45
gen aged45 = 1
replace aged45 = 0 if age_filing>=1 & age_filing<45
tab aged45

* yrs0609:	0 "1997-2005"	1 "2006-2009"
cap drop 
gen yrs0609=0
replace yrs0609=1 if startpay_year>=2006 & startpay_year<=2009
tab yrs0609

*** Variable descriptions - other identifiers

* primary disability: individual dummy variables identifying them
* "mental"
* "musc"
* "neoplasms"
* "circulatory"
* "other"

*** RKD with covariates for 0/1 identifiers

local vlist "male black aged45 yrs0609 mental musc neoplasms circulatory other"
foreach v of local vlist  {	

forvalues i=0/1	{
	
preserve
keep if `v'==`i'

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"
	
replace categ = "`v'_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace

restore

}

}
	
*** Mortality in first 2 years

foreach i in 12 34	{
 
use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* DDS restriction
keep if hearings==0

* time on DI - years 1&2, years 3&4
cap drop d_year12
gen d_year12 = (d_year1 + d_year2) / 2
cap drop d_year34
gen d_year34 = (d_year3 + d_year4) / 2

* Fuzzy RK with covariates 
capture noisily rdrobust d_year`i' ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"

regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year`i'

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

replace categ = "mortality_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}

*** Adding column and group identifiers

use "${output}/tables/`tbl'.dta", clear

cap drop N
cap drop if var=="RD_Estimate"
cap drop if coef==.
cap drop _id

cap gen column=3
replace column=4 if var=="mean" 

gen obsnum = _n
gsort column -obsnum

drop obsnum

order categ
order column

save "${output}/tables/`tbl'.dta", replace

}


**************************************************************************************
**************************************************************************************
********************************** APPENDIX TABLES ***********************************
**************************************************************************************
**************************************************************************************

**************************************************************************************
************************************* TABLE A2 ***************************************
**************************************************************************************
*** Summary statistics

local tbl tablea2
local sample main

*** Summary statistics for each sample

local bplist lower fm upper
foreach bp of local bplist  {
	
use "${dta}/`bp'_`sample'_sample.dta", clear

*Annualized payments
gen pia_annual = pia*12
gen fmax_annual = fmax*12
*Dependents identifer 
gen depend=(family_ben==1 & (diff_start==-1 | diff_start==0))

keep if imeabs<=1200

#delimit cr
summ age_filing male black 
	pia pia_annual depend fmax fmax_annual 
	hearings musc mental neoplasms circulatory other 
	d_year1 d_year2 d_year3 d_year4;
#delimit cr

count
  
}  

*** Summary statistics overall
	
use "${dta}\master_file2.dta", clear

*Annualized payments
gen pia_annual = pia*12
gen fmax_annual = fmax*12
*Dependents identifer 
gen depend=(family_ben==1 & (diff_start==-1 | diff_start==0))

#delimit cr
summ age_filing male black 
	pia pia_annual depend fmax fmax_annual 
	hearings musc mental neoplasms circulatory other 
	d_year1 d_year2 d_year3 d_year4;
#delimit cr

count


**************************************************************************************
************************************* TABLE B1 ***************************************
**************************************************************************************
*** Testing for kinks in the densities at the bend points

local tbl tableb1

local sample main

*** Creating a data set to store the regression results 
preserve
clear
set obs 1
gen kink_pval = .
save "${output}/tables/`tbl'.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {
	
use "${dta}/`bp'_`sample'_sample.dta", clear
gen counter=1
gen mid=(ceil(ime)) - 0.5
collapse (sum) counter, by(mid)
rename mid ime

cap gen bendpoint=""
replace bendpoint="`bp'" if bendpoint==""

if bendpoint=="lower" local restrict = 194.48
if bendpoint=="fm"    local restrict = 599.573
if bendpoint=="upper" local restrict = 588.943

if bendpoint=="lower" local binsize = 194.48 /8
if bendpoint=="fm"    local binsize = 599.573/24
if bendpoint=="upper" local binsize = 588.943/24

expand counter

keep ime

*** Restrictions
cap drop imeabs
gen imeabs=abs(ime)
keep if imeabs <= `restrict'

count
local sample_size = r(N)

local max_order = 3

mat mat_pvals = J(`max_order',4,.)

qui tempvar d x
	
qui gen `d' = ime >= 0
	
qui gen `x' = floor(ime/(`binsize'))*(`binsize')+0.5*(`binsize')
	
qui bys `x': gen frac = (_N/`sample_size')

qui bys `x': gen dens_est = frac/(`binsize')
	
qui bys `x': keep if _n==1
	
qui gen double weight = 1/frac

** Note: specify the set of polynomial order candidates
** for fitting the density and testing a kink

forval fit_order = 1/`max_order' {

	di ""
	di "Density estimation with order `fit_order'" 

	preserve
	
	local regvar_full = ""	
	
	forval k=1(1)`fit_order' {
		
		qui cap drop x`k'
		qui cap drop dx`k'
		
		qui gen x`k'=`x'^(`k')
		qui gen dx`k'=`d'*`x'^(`k')

		local regvar_full = "`regvar_full'" + "x`k' " + "dx`k' "
		
						}							
	
	** regvar_full includes the d term
	qui gen d = `d'
	local regvar_full = "`regvar_full'" + "d"
	
	** regvar_r excludes the dx1 as well as the d term
	local regvar_r = subinstr("`regvar_full'", "dx1 ", "", .)
	
	reg frac `regvar_full' [aweight=weight]
	qui predict frac_hat_full
	qui predict res_full, res
	local gof_full_df = e(df_r)-1
	
	qui gen frac_hat_full_left = frac_hat_full if `x'<0
	qui replace frac_hat_full_left = frac_hat_full - _b[d] if `x'==0
	
	reg frac `regvar_r' [aweight=weight]
	qui predict frac_hat_r
	qui predict res_r, res
	
	qui gen dev_full = res_full^2*weight
	qui gen dev_r = res_r^2*weight
	
	qui total dev_full
	local pchi2_full = _b[dev_full] * `sample_size'
	
	qui total dev_r
	local pchi2_r = _b[dev_r] * `sample_size'
	
	** p_values
	local kink_pval = 1-chi2(1,`pchi2_r'-`pchi2_full')
	local gof_full_pval = 1-chi2(`gof_full_df',`pchi2_full')
	
	mat mat_pvals[`fit_order',1] = `fit_order'
	mat mat_pvals[`fit_order',2] = `kink_pval'
	mat mat_pvals[`fit_order',3] = `gof_full_pval'	
	mat mat_pvals[`fit_order',4] = `pchi2_full'	
	
	restore		
	}

clear 

svmat mat_pvals, n("var")
rename var1 order
rename var2 kink_pval
rename var3 gof_full_pval
rename var4 pchi2_full

gen pchi2_full_last = pchi2_full[_n-1]
gen pchi2_full_diff = pchi2_full_last - pchi2_full
gen gof_full_rel_pval = 1 - chi2(2,pchi2_full_diff)

keep order kink_pval

cap gen bendpoint=""
replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

order column
order bendpoint

append using "${output}/tables/`tbl'.dta"

drop if kink_pval==.
sort column order

save "${output}/tables/`tbl'.dta", replace

}


**************************************************************************************
************************************* TABLE B2 ***************************************
**************************************************************************************
*** Estimated kinks in beneficiary characteristics 

local tbl tableb2

local sample main

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
gen categ=""
save "${output}/tables/`tbl'.dta", replace
restore

*** Generating the regression results - main covariates

local bplist lower fm upper
foreach bp of local bplist  {

local vlist "age_filing male black hearings mental musc"
foreach v of local vlist  {	

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* RK with covariates 
capture noisily rdrobust `v' ime, c(0) deriv(1) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) 

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc)
  local se_rob = e(se_tau_rb)

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef="`v'" 

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

cap drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""
replace categ = "`v'" if categ==""

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}

}

*** Generating the regression results - SSI

local sample ssi

local bplist lower fm upper
foreach bp of local bplist  {

local vlist "ssi"
foreach v of local vlist  {	

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/
	
* RK with covariates 
capture noisily rdrobust `v' ime, c(0) deriv(1) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) 

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) 
  local se_rob = e(se_tau_rb)

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef="`v'" 

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

cap drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""
replace categ = "`v'" if categ==""

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}

}

*** Adding column and group identifiers

use "${output}/tables/`tbl'.dta", clear

cap gen rorder=.
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

cap gen column=.
replace column=1 if categ=="age_filing"
replace column=2 if categ=="male"
replace column=3 if categ=="black"
replace column=4 if categ=="hearings"
replace column=5 if categ=="mental"
replace column=6 if categ=="musc"
replace column=7 if categ=="ssi"

cap gen panel = .
replace panel = 1 if bendpoint=="lower"
replace panel = 2 if bendpoint=="fm"
replace panel = 3 if bendpoint=="upper"

sort panel column rorder
drop _id
drop rorder

order var
order categ
order column
order panel
order bendpoint

save "${output}/tables/`tbl'.dta", replace


**************************************************************************************
************************************* TABLE B3 ****************************************
**************************************************************************************
*** Placebo estimates

local tbl tableb3

local cov_controls = "age_filing male black hearings musc mental"

**************************************************************************************
*** Top panel: DI beneficiaries without dependents (Family bend point only)

local sample nondepend

local bp fm

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'_`sample'.dta", replace
restore

*** Sharp RK with covariates - no first stage in DI payments

use "${dta}/fm_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

capture noisily rdrobust d_year14 ime, c(0) deriv(1) scalepar(22522.5) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - already converted to p.p. change per $1000/year by scalepar
  local b_rob = e(tau_bc)
  local se_rob = e(se_tau_rb)

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

*** column 2 (no columns 1 and 3)
regsave using "${output}/tables/`tbl'_`sample'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14 

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'_`sample'.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

sort column rorder
drop _id
drop rorder

order var
order column
order bendpoint

save "${output}/tables/`tbl'_`sample'.dta", replace


**************************************************************************************
*** Middle panel: DI beneficiaries' predicted mortality based on predetermined covariates

local sample predict

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'_`sample'.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_main_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

*** Create covariate index
xtile age_q = age_filing, nq(10)
reg d_year14 i.age_q male black hearings musc mental
predict covind_full

*** Fuzzy RK with covariates 
capture noisily rdrobust covind_full ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

regsave using "${output}/tables/`tbl'_`sample'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=covind_full

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'_`sample'.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

sort column rorder
drop _id
drop rorder

order var
order column
order bendpoint

save "${output}/tables/`tbl'_`sample'.dta", replace

}


**************************************************************************************
*** Bottom panel: Non-beneficiaries

local sample nonbenef

local cov_controls2 = "age_cwhs male black" /*There are a more limited set of controls available in these data*/

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'_`sample'.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

cap gen bendpoint=""
replace bendpoint="`bp'" if bendpoint==""

if bendpoint=="lower" local scale = 14367.8
if bendpoint=="fm"    local scale = 22522.5
if bendpoint=="upper" local scale = 49019.6

*** Sharp RK with covariates 
capture noisily rdrobust mort ime, c(0) deriv(1) scalepar(`scale') p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls2')

  *main results - already converted to p.p. change per $1000/year by scalepar
  local b_rob = e(tau_bc)
  local se_rob = e(se_tau_rb)

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

regsave using "${output}/tables/`tbl'_`sample'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=mort

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'_`sample'.dta"

sort _id

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

sort column rorder
drop _id
drop rorder

order var
order column
order bendpoint

save "${output}/tables/`tbl'_`sample'.dta", replace

}


**************************************************************************************
************************************* TABLE B4 ***************************************
**************************************************************************************
*** Table B4 Heterogeneity in the mortality effects at the upper bend point 

local bp upper

local sample main

**************************************************************************************
*** Full sample

local tbl tableb4_all

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef = .
gen categ = ""
save "${output}/tables/`tbl'.dta", replace
restore

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

*** Variable descriptions - 0/1 identifiers

* hearings: 0 "DDS"	 		1 "hearings"

* male: 	0 "females"		1 "males"

* black: 	0 "nonblack"	1 "black"

* aged45:	0 "aged<45"		1 "aged 45+"
cap drop aged45
gen aged45 = 1
replace aged45 = 0 if age_filing>=1 & age_filing<45
tab aged45

* yrs0609:	0 "1997-2005"	1 "2006-2009"
cap drop 
gen yrs0609=0
replace yrs0609=1 if startpay_year>=2006 & startpay_year<=2009
tab yrs0609

*** Variable descriptions - other identifiers

* primary disability: individual dummy variables identifying them
* "mental"
* "musc"
* "neoplasms"
* "circulatory"
* "other"

*** RKD with covariates for 0/1 identifiers

local vlist "hearings male black aged45 yrs0609 mental musc neoplasms circulatory other"
foreach v of local vlist  {	

forvalues i=0/1	{
	
preserve
keep if `v'==`i'

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"
	
replace categ = "`v'_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace

restore

}

}

*** Mortality in first 2 years

foreach i in 12 34	{
 
use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* time on DI - years 1&2, years 3&4
cap drop d_year12
gen d_year12 = (d_year1 + d_year2) / 2
cap drop d_year34
gen d_year34 = (d_year3 + d_year4) / 2

* Fuzzy RK with covariates 
capture noisily rdrobust d_year`i' ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"

regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year`i'

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

replace categ = "mortality_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}

*** Adding column and group identifiers

use "${output}/tables/`tbl'.dta", clear

cap drop N
cap drop if var=="RD_Estimate"
cap drop if coef==.
cap drop _id

cap gen column=1
replace column=2 if var=="mean" 

gen obsnum = _n
gsort column -obsnum

drop obsnum

order categ
order column

save "${output}/tables/`tbl'.dta", replace

**************************************************************************************
*** DDS only

local tbl tableb4_dds

local cov_controls = "age_filing male black musc mental"

*** Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef = .
gen categ = ""
save "${output}/tables/`tbl'.dta", replace
restore

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* DDS restriction
keep if hearings==0

*** Variable descriptions - 0/1 identifiers

* hearings: 0 "DDS"	 		1 "hearings"

* male: 	0 "females"		1 "males"

* black: 	0 "nonblack"	1 "black"

* aged45:	0 "aged<45"		1 "aged 45+"
cap drop aged45
gen aged45 = 1
replace aged45 = 0 if age_filing>=1 & age_filing<45
tab aged45

* yrs0609:	0 "1997-2005"	1 "2006-2009"
cap drop 
gen yrs0609=0
replace yrs0609=1 if startpay_year>=2006 & startpay_year<=2009
tab yrs0609

*** Variable descriptions - other identifiers

* primary disability: individual dummy variables identifying them
* "mental"
* "musc"
* "neoplasms"
* "circulatory"
* "other"

*** RKD with covariates for 0/1 identifiers

local vlist "male black aged45 yrs0609 mental musc neoplasms circulatory other"
foreach v of local vlist  {	

forvalues i=0/1	{
	
preserve
keep if `v'==`i'

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"
	
replace categ = "`v'_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace

restore

}

}
	
*** Mortality in first 2 years

foreach i in 12 34	{
 
use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* DDS restriction
keep if hearings==0

* time on DI - years 1&2, years 3&4
cap drop d_year12
gen d_year12 = (d_year1 + d_year2) / 2
cap drop d_year34
gen d_year34 = (d_year3 + d_year4) / 2

* Fuzzy RK with covariates 
capture noisily rdrobust d_year`i' ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"

regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(m_rob, `b_rob', `se_rob') append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year`i'

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

replace categ = "mortality_`i'" if categ==""	

order var
order categ

save "${output}/tables/`tbl'.dta", replace
}


*** Adding column and group identifiers

use "${output}/tables/`tbl'.dta", clear

cap drop N
cap drop if var=="RD_Estimate"
cap drop if coef==.
cap drop _id

cap gen column=3
replace column=4 if var=="mean" 

gen obsnum = _n
gsort column -obsnum

drop obsnum

order categ
order column

save "${output}/tables/`tbl'.dta", replace


**************************************************************************************
************************************* TABLE B5 ***************************************
**************************************************************************************
*** Effect of DI Benefits on mortality over longer follow-up periods 

local tbl tableb5

local sample main

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'.dta", replace
restore

*** Generating the regression results

local bplist lower fm upper
foreach bp of local bplist  {

forvalues i=4/10	{
	
use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* Fuzzy RK with covariates 
capture noisily rdrobust d_year1`i' ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *first-stage results
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 

keep if ime>=-10 & ime<0
collapse (mean) coef=d_year1`i' 

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen mort_yrs = .
replace mort_yrs = `i' if mort_yrs == .

cap gen column = .
replace column = `i' - 3 if column==.

cap gen rorder=.
replace rorder=0 if var=="f_rob"
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

cap gen panel = .
replace panel = 1 if bendpoint=="lower"
replace panel = 2 if bendpoint=="fm"
replace panel = 3 if bendpoint=="upper"

sort panel column rorder
drop _id
drop rorder

order var
order mort_yrs
order column
order panel
order bendpoint

save "${output}/tables/`tbl'.dta", replace

}

}

**************************************************************************************
************************************* TABLE B6 ***************************************
**************************************************************************************
*** Effect of DI payments on earnings

local tbl tableb6

local sample main

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'.dta", replace
restore

*** Generating the regression results - MSE-optimal bandwidths

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* Fuzzy RK without covariates 
capture noisily rdrobust e1234 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3)  

  *main results - converted to p.p. change per $1000/year
  local b_conv = e(tau_cl) * (1/12) *(1000)* 100
  local se_conv = e(se_tau_cl) * (1/12) * 1000  * 100
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Conv coefficient: " "`b_conv'"
  di "Conv std error: "   "`se_conv'"
  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Conventional 
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_conv, `b_conv', `se_conv', bw_conv, `bw', miss, obs_conv, `obs', miss) append

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid /// 
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append


* Fuzzy RK with covariates 
capture noisily rdrobust e1234 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Generating the regression results - bandwidths from Gelber, Moore & Strand (2017)
use "${dta}/`bp'_`sample'_sample.dta", clear
	
cap gen bendpoint=""
replace bendpoint="`bp'" if bendpoint==""

if bendpoint=="lower" local bwidth = 600
if bendpoint=="fm"    local bwdith = 900
if bendpoint=="upper" local bwidth = 1500

*** Fuzzy RK  
capture noisily rdrobust e1234 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) h(`bwidth') kernel(uni) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_conv = e(tau_cl) * (1/12) *(1000)* 100
  local se_conv = e(se_tau_cl) * (1/12) * 1000  * 100
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Conv coefficient: " "`b_conv'"
  di "Conv std error: "   "`se_conv'"
  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Conventional 
regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(m_conv, `b_conv', `se_conv', bw_conv, `bw', miss, obs_conv, `obs', miss) append

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid /// 
	addvar(m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Fuzzy RK with covariates 
capture noisily rdrobust e1234 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) h(`bwidth') kernel(uni) vce(nn 3) covs(`cov_controls')

  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

* Robust
regsave using "${output}/tables/`tbl'.dta", autoid  ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Organizing the results 
use "${output}/tables/`tbl'.dta", clear

drop N
drop if var=="RD_Estimate"
drop if coef==.

cap gen column = _id
replace column = _id if column==.
drop _id

replace bendpoint="`bp'" if bendpoint==""

cap gen rorder=.
replace rorder=1 if var=="m_conv"
replace rorder=2 if var=="bw_conv"
replace rorder=3 if var=="obs_conv"
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"

cap gen panel = .
replace panel = 1 if bendpoint=="lower"
replace panel = 2 if bendpoint=="fm"
replace panel = 3 if bendpoint=="upper"

sort panel column rorder
drop rorder

order var
order column
order panel
order bendpoint

save "${output}/tables/`tbl'.dta", replace

}

drop if column>3 & bendpoint=="fm"

save "${output}/tables/`tbl'.dta", replace


**************************************************************************************
************************************* TABLE B7 ****************************************
**************************************************************************************
*** Combined effect of DI payments plus earnings on mortality rates

local tbl tableb7

local sample main

local cov_controls = "age_filing male black hearings musc mental"

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

*** Add the earnings to DI payments
gen double fmax_e = fmax + e1234
replace fmax_e = fmax if fmax_e==.

summ fmax_e

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax_e) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *first-stage results
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

regsave using "${output}/tables/`tbl'.dta", autoid ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=0 if var=="f_rob"
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

cap gen panel = .
replace panel = 1 if bendpoint=="lower"
replace panel = 2 if bendpoint=="fm"
replace panel = 3 if bendpoint=="upper"

sort panel column rorder
drop _id
drop rorder

order var
order column
order panel
order bendpoint

save "${output}/tables/`tbl'.dta", replace

}

**************************************************************************************
************************************* TABLE B8 ****************************************
**************************************************************************************
*** Estimates using DI/SSI beneficiaries after the DI waiting period

local tbl tableb8

local sample ssi

local cov_controls = "age_filing male black hearings musc mental"

**************************************************************************************
*** Top panel: DI beneficiaries who received SSI

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'_`sample'_only.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

* Restrict to only SSI
keep if ssi==1

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *first-stage results
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

regsave using "${output}/tables/`tbl'_`sample'_only.dta", autoid ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'_`sample'_only.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=0 if var=="f_rob"
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

sort column rorder
drop _id
drop rorder

order var
order column
order bendpoint

save "${output}/tables/`tbl'_`sample'_only.dta", replace

}


**************************************************************************************
*** Bottom panel: All DI beneficiaries – Main sample + DI beneficiaries who received SSI

*** Creating a data set to store the regression results 

preserve
clear
set obs 1
gen coef=.
gen bendpoint=""
save "${output}/tables/`tbl'_`sample'_all.dta", replace
restore

local bplist lower fm upper
foreach bp of local bplist  {

use "${dta}/`bp'_`sample'_sample.dta", clear
keep if imeabs<=1200 /*common starting data*/

*** Fuzzy RK with covariates 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

  *first-stage results
  local f_b_rob = e(tau_T_bc)
  local f_se_rob = e(se_tau_T_cl)
  
  *main results - converted to p.p. change per $1000/year
  local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
  local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

  di "Robust coefficient: " "`b_rob'"
  di "Robust std error: "   "`se_rob'"
  
  *extra information - bandwidth and observations
  local bw = e(h_l)
  local obs = e(N_h_l) + e(N_h_r)
  
  di "Bandwidth: "    "`bw'"
  di "Observations: " "`obs'"

regsave using "${output}/tables/`tbl'_`sample'_all.dta", autoid ///
	addvar(f_rob, `f_b_rob', `f_se_rob', m_rob, `b_rob', `se_rob', bw_rob, `bw', miss, obs_rob, `obs', miss) append

*** Mean within $10 below the bend point 
keep if ime>=-10 & ime<0
collapse (mean) coef=d_year14

gen var="mean"
keep coef var

append using "${output}/tables/`tbl'_`sample'_all.dta"

drop N
drop if var=="RD_Estimate"
drop if coef==.

replace bendpoint="`bp'" if bendpoint==""

cap gen column = .
replace column = 1 if bendpoint=="lower"
replace column = 2 if bendpoint=="fm"
replace column = 3 if bendpoint=="upper"

cap gen rorder=.
replace rorder=0 if var=="f_rob"
replace rorder=1 if var=="m_rob"
replace rorder=2 if var=="bw_rob"
replace rorder=3 if var=="obs_rob"
replace rorder=4 if var=="mean"

sort column rorder
drop _id
drop rorder

order var
order column
order bendpoint

save "${output}/tables/`tbl'_`sample'_all.dta", replace

}

clear all

log close